Summary

This script contains code to reproduce the analyses and figures from Hinchliffe C et al 2020, “Larval fish distribution in a changing western boundary current” (Link).

There are 3 major analyses within the manuscript. The first two involve analysis of the seasonal and latitudinal variation in larval fish abundance and taxa richness using generalised additive mixed-models (GAMMs) and can be found in section 2) and 3) in this script. The third major analyses involves investigation of the community composition of larval fish across the study area, first in relation to Sea Surface Temperature (SST), followed by Latitude (which appears in the supplementary material of the manuscript), and can be found in section 4) of this script. The figures associated with each analysis are at the end of each section, with there figure number as they appear in the manuscript and figure caption.

1) Data manipulation and cleaning

This section of the script cleans data and creates variables required for the subsequent analyses.

load required packages for all analyses and figure building.

library(tidyverse)
library(mgcv)
library(vegan)
library(car)
library(lubridate)
library(pracma)
#install_github('gavinsimpson/gratia@v0.2-3')
library(gratia)
library(lctools)
library(ggplot2)
library(gridGraphics)
library(cowplot)
library(ggpubr)
library(grid)
library(ggplotify)
library(mvabund)
#install ecoCopula
#devtools::install_github("gordy2x/ecoCopula",force=TRUE)
library(ecoCopula)
library(quantreg)

bring in data

nimo.data <- read.csv("data/all_data.csv", header=TRUE)
#head(nimo.data)

Make bathymetry absolute values in km units

nimo.data$Bathy <- (abs(nimo.data$Bathy)/1000)   
#summary(nimo.data$Bathy)

Create date variables

nimo.data$DateCaught <- as.Date(nimo.data$Date, "%d/%m/%Y")
nimo.data$Month <- month(nimo.data$Date)
#summary(nimo.data$Month)
nimo.data$Year <- year(nimo.data$DateCaught)
nimo.data$Day <- day(nimo.data$DateCaught)
#summary(nimo.data$Day)

nimo.data$dayofyear <- yday(nimo.data$DateCaught)
#summary(nimo.data$dayofyear)
#hist(nimo.data$dayofyear)

Create season variables

nimo.data$Month <- as.factor(nimo.data$Month) #month as factor

#make season variable
nimo.data$Season <- NA_character_ 

nimo.data$Season[nimo.data$Month %in% c(1,2,12)] <- "summer"
nimo.data$Season[nimo.data$Month %in% 3:5] <- "autumn"
nimo.data$Season[nimo.data$Month %in% 6:8] <- "winter"
nimo.data$Season[nimo.data$Month %in% 9:11] <- "spring"

nimo.data$Season <- as.factor(nimo.data$Season)
nimo.data$Month <- as.numeric(nimo.data$Month) 

Create time period variables

nimo.data$Decade[nimo.data$Year  <= 1999] <- "1990s"
nimo.data$Decade[nimo.data$Year  <= 2009 & nimo.data$Year > 1999] <- "2000s"
nimo.data$Decade[nimo.data$Year  <= 2020 & nimo.data$Year > 2009] <- "2010s"

nimo.data$Period[nimo.data$Year  <= 1998] <- "pre1998"
nimo.data$Period[nimo.data$Year  >= 1999] <- "post1998"

nimo.data$Decade <- as.factor(nimo.data$Decade)
#summary(nimo.data$Decade)

nimo.data$Period <- as.factor(nimo.data$Period)
#summary(nimo.data$Period)

Create latitudinal region variable

nimo.data$Region[nimo.data$Latitude  <= -40] <- "South"

nimo.data$Region[nimo.data$Latitude  <= -35 & nimo.data$Latitude > -40] <- "Mid-South"

nimo.data$Region[nimo.data$Latitude  <= -30 & nimo.data$Latitude > -35] <- "Mid-North"

nimo.data$Region[nimo.data$Latitude  > -30] <- "North"

nimo.data$Region <- as.factor(nimo.data$Region)

nimo.data$Region <- as.factor(nimo.data$Region)
nimo.data$Region <- factor(nimo.data$Region, ordered = TRUE, levels = c("North", "Mid-North", "Mid-South","South"))
#summary(nimo.data$Region) 

Create mid-point standardised sampling depth variable

nimo.data$stand_depth <- as.character(nimo.data$Gear_depth_m)



nimo.data$stand_depth[nimo.data$stand_depth  == "0-25"] <- "25" #double mid point

nimo.data$stand_depth[nimo.data$stand_depth  == "0-30"] <- "30"

nimo.data$stand_depth[nimo.data$stand_depth  == "0-35"] <- "35"

nimo.data$stand_depth[nimo.data$stand_depth  == "0-40"] <- "40"

nimo.data$stand_depth[nimo.data$stand_depth  == "0-50"] <- "50"

nimo.data$stand_depth[nimo.data$stand_depth  == "0-80"] <- "80"

nimo.data$stand_depth[nimo.data$stand_depth  == "25-50"] <- "75"

nimo.data$stand_depth[nimo.data$stand_depth  == "50-75"] <- "125"

nimo.data$stand_depth[nimo.data$stand_depth  == "75-100"] <- "175"

nimo.data$stand_depth <- as.factor(nimo.data$stand_depth)

nimo.data$stand_depth <- as.numeric(nimo.data$stand_depth)

nimo.data$stand_depth <- (nimo.data$stand_depth/2)

Create sampling depth factors

nimo.data$depth_sect[nimo.data$stand_depth  >= 0 & nimo.data$stand_depth < 11] <- "a) 0-10"

nimo.data$depth_sect[nimo.data$stand_depth  >= 11 & nimo.data$stand_depth < 50] <- "b) 11-50"

nimo.data$depth_sect[nimo.data$stand_depth  >= 51 & nimo.data$stand_depth < 100] <- "c) 51-100"

nimo.data$depth_sect[nimo.data$stand_depth  >= 101 & nimo.data$stand_depth < 150] <- "d) 101-150"

nimo.data$depth_sect[nimo.data$stand_depth  >= 151] <- "e) 151+"

nimo.data$depth_sect <- as.factor(nimo.data$depth_sect)

Creating diversity variables

larval.spp<-nimo.data[1:3193,22:241]

abundance <- rowSums(larval.spp)
richness <- specnumber(larval.spp) # richness
shannon <- diversity(larval.spp) #shannon weaver index from vegan package
J <- shannon/log(specnumber(larval.spp)) # Evenness
shannon_effective <- exp(diversity(larval.spp))


nimo.data$shannon <- shannon
nimo.data$evenness <- J
nimo.data$richness <- richness
nimo.data$abundance <- abundance
nimo.data$shannon_effective <- shannon_effective

Separate Australian west coast and east coast

Samples into separate data frames

eastcoast <- subset(nimo.data, nimo.data$Longitude >= 130)
westcoast <- subset(nimo.data, nimo.data$Longitude <= 130)

Check SST and Hadley Centre Global Sea Ice and Sea Surface Temperature (HadISST) correlation on east coast

eastcoast$SST <- as.numeric(eastcoast$SST)
eastcoast$HadISST <- as.numeric(eastcoast$HadISST)

cor.test(eastcoast$SST,eastcoast$HadISST, method = c("pearson", "kendall", "spearman"))
## 
##  Pearson's product-moment correlation
## 
## data:  eastcoast$SST and eastcoast$HadISST
## t = 108.89, df = 1443, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9383103 0.9495273
## sample estimates:
##       cor 
## 0.9441919
#plot(x = eastcoast$SST,y = eastcoast$HadISST)

Separate inshore and offshore

Sites on the east coast based on position on continental shelf (boundary set at 200m bathymetry).

inshore <- subset(eastcoast, eastcoast$Bathy <= .200)
offshore <- subset(eastcoast, eastcoast$Bathy >= .200)

2) Abundance GAMM

The following section contains code for the analysis of larval fish abundance using generalised additive mixed models (GAMM).

Testing linear correlation between predictors using Pearson correlation

preds <- data.frame(eastcoast$Latitude, eastcoast$HadISST, eastcoast$stand_depth, eastcoast$Bathy, eastcoast$dists_km, eastcoast$Temperature_C, eastcoast$SST)
cor(preds)
##                         eastcoast.Latitude eastcoast.HadISST eastcoast.stand_depth eastcoast.Bathy eastcoast.dists_km eastcoast.Temperature_C eastcoast.SST
## eastcoast.Latitude               1.0000000         0.7939551           -0.20202890              NA         0.29394633                      NA            NA
## eastcoast.HadISST                0.7939551         1.0000000           -0.16558965              NA         0.27806479                      NA            NA
## eastcoast.stand_depth           -0.2020289        -0.1655896            1.00000000              NA         0.00426594                      NA            NA
## eastcoast.Bathy                         NA                NA                    NA               1                 NA                      NA            NA
## eastcoast.dists_km               0.2939463         0.2780648            0.00426594              NA         1.00000000                      NA            NA
## eastcoast.Temperature_C                 NA                NA                    NA              NA                 NA                       1            NA
## eastcoast.SST                           NA                NA                    NA              NA                 NA                      NA             1

Model selection

# step 1) Full model
model1 <-  gam(abundance ~  Season + Period +  exp(-Bathy)  +  s(stand_depth, bs =  "cs", k = 5) + s(dists_km, bs =  "cs", k = 5) + s(HadISST, bs =  "cs", k = 5)  + s(Latitude, bs =  "fs", xt = "cs", k = 5, by = Season, id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re"), data = eastcoast, family = nb(theta = NULL, link = "log"), method = "ML", fit=TRUE,  Select=TRUE)


concurvity(model1, full=T) #checking concurvity
##          para s(stand_depth) s(dists_km) s(HadISST) s(Latitude):Seasonautumn s(Latitude):Seasonspring s(Latitude):Seasonsummer s(Latitude):Seasonwinter s(Project_name)
## worst       1      0.3174216   0.8765001  0.9479951                0.9191817                0.9259444                0.9537419                0.9044561       1.0000000
## observed    1      0.1266766   0.7673233  0.9215439                0.8004886                0.8111269                0.9347195                0.7972925       0.9589919
## estimate    1      0.1626038   0.6044313  0.8125066                0.7369388                0.7975847                0.8583453                0.8191350       0.8065563
# step 2)  s(dists_km, bs =  "cs", k = 5) + s(HadISST, bs =  "cs", k = 5)  removed due to concurvity
model2 <-  gam(abundance ~  Season + Period +  exp(-Bathy)  +  s(stand_depth, bs =  "cs", k = 5)  + s(Latitude, bs =  "fs", xt = "cs", k = 5, by = Season, id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re"), data = eastcoast, family = nb(theta = NULL, link = "log"), method = "ML", fit=TRUE,  Select=TRUE)


# Step 3) Period removed based on AIC score, checked with chi-sq test
model2i <-  gam(abundance ~  Season  +  exp(-Bathy)  +  s(stand_depth, bs =  "cs", k = 5)  + s(Latitude, bs =  "fs", xt = "cs", k = 5, by = Season, id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re"), data = eastcoast, family = nb(theta = NULL, link = "log"), method = "ML", fit=TRUE,  Select=TRUE)

AIC(model2,model2i) ## AIC reduced with period removed
##               df      AIC
## model2  30.69146 32233.64
## model2i 29.84881 32231.60
anova(model2,model2i, test = "Chisq") #non sig therefore can remove period (adds nothing to model)
## Analysis of Deviance Table
## 
## Model 1: abundance ~ Season + Period + exp(-Bathy) + s(stand_depth, bs = "cs", 
##     k = 5) + s(Latitude, bs = "fs", xt = "cs", k = 5, by = Season, 
##     id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re")
## Model 2: abundance ~ Season + exp(-Bathy) + s(stand_depth, bs = "cs", 
##     k = 5) + s(Latitude, bs = "fs", xt = "cs", k = 5, by = Season, 
##     id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re")
##   Resid. Df Resid. Dev       Df Deviance Pr(>Chi)
## 1    2968.1     3580.7                           
## 2    2969.0     3580.7 -0.93486 0.021107
# Step 4) linear Season term removed, but increases AIC and is sig dif in chi-sq test there retain Season in model
model3 <-  gam(abundance ~  exp(-Bathy)  +  s(stand_depth, bs =  "cs", k = 5)  + s(Latitude, bs =  "fs", xt = "cs", k = 5, by = Season, id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re"), data = eastcoast, family = nb(theta = NULL, link = "log"), method = "ML", fit=TRUE,  Select=TRUE)


AIC(model3,model2i) ## AIC increased when season removed. ## model2i best model here. 
##               df      AIC
## model3  29.40232 32256.84
## model2i 29.84881 32231.60
anova(model3,model2i, test = "Chisq")  #sig dif therefore dont remove season
## Analysis of Deviance Table
## 
## Model 1: abundance ~ exp(-Bathy) + s(stand_depth, bs = "cs", k = 5) + 
##     s(Latitude, bs = "fs", xt = "cs", k = 5, by = Season, id = 1) + 
##     offset(log(Volume_m3)) + s(Project_name, bs = "re")
## Model 2: abundance ~ Season + exp(-Bathy) + s(stand_depth, bs = "cs", 
##     k = 5) + s(Latitude, bs = "fs", xt = "cs", k = 5, by = Season, 
##     id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re")
##   Resid. Df Resid. Dev      Df Deviance Pr(>Chi)  
## 1    2969.4     3582.7                            
## 2    2969.0     3580.7 0.41619   1.9727  0.05599 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## therefore best model for abundance appears to be model2i, removing Period based on AIC and s(dists_km, bs =  "cs", k = 5) + s(HadISST, bs =  "cs", k = 5) based on concurvity.

#Run with method = REML for estimates
model2i <-  gam(abundance ~  Season  +  exp(-Bathy)  +  s(stand_depth, bs =  "cs", k = 5)  + s(Latitude, bs =  "fs", xt = "cs", k = 5, by = Season, id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re"), data = eastcoast, family = nb(theta = NULL, link = "log"), method = "REML", fit=TRUE,  Select=TRUE)

#Best model summary
summary(model2i) 
## 
## Family: Negative Binomial(0.586) 
## Link function: log 
## 
## Formula:
## abundance ~ Season + exp(-Bathy) + s(stand_depth, bs = "cs", 
##     k = 5) + s(Latitude, bs = "fs", xt = "cs", k = 5, by = Season, 
##     id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re")
## 
## Parametric coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.93971    0.30561  -6.347 2.19e-10 ***
## Seasonspring -0.55675    0.11191  -4.975 6.52e-07 ***
## Seasonsummer  0.11398    0.06852   1.664  0.09621 .  
## Seasonwinter -0.31556    0.10983  -2.873  0.00406 ** 
## exp(-Bathy)   0.96264    0.13170   7.309 2.69e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                            edf Ref.df Chi.sq  p-value    
## s(stand_depth)           3.866      4  582.6  < 2e-16 ***
## s(Latitude):Seasonautumn 2.377      4  303.0  < 2e-16 ***
## s(Latitude):Seasonspring 1.967      4  442.4 3.92e-08 ***
## s(Latitude):Seasonsummer 2.469      4  817.2 5.58e-10 ***
## s(Latitude):Seasonwinter 2.064      4  273.9 1.96e-08 ***
## s(Project_name)          9.487     10  230.2  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  -0.079   Deviance explained = 30.2%
## -REML =  16152  Scale est. = 1         n = 3000
#par(mfcol=c(3,2))
#plot(model2i, shade = TRUE)


#checking concurvity
#concurvity(model2i,full=TRUE)

#checking assumptions
appraise(model2i)

#over disperison test
E1 <- resid(model2i, type = "pearson")   #ok - less than 2
N <- nrow(eastcoast)
p <- length(coef(model2i))
Overdispersion <- sum(E1^2) / (N - p)
#paste('Overdispersion: ',Overdispersion) #"Overdispersion:  1.58845886919077"

Checking spatial auto correlation of Abundance GAMM with Moran’s I If Moran’s I is ~0 then there is no spatial autocorrelation.

res <- residuals(model2i)
#res


Coords <-cbind(eastcoast$Latitude, eastcoast$Longitude)

bw <- 6

mI <-moransI(Coords,bw,res)

moran.table <-matrix(data=NA,nrow=1,ncol=6)

col.names <-c("Moran's I", "Expected I", "Z resampling", "P-value resampling","Z randomization", "P-value randomization")

colnames(moran.table) <- col.names

moran.table[1,1] <- mI$Morans.I

moran.table[1,2] <- mI$Expected.I

moran.table[1,3] <- mI$z.resampling

moran.table[1,4] <- mI$p.value.resampling

moran.table[1,5] <- mI$z.randomization

moran.table[1,6] <- mI$p.value.randomization

print(moran.table)

Plotting output

Plotting effects of the larval fish abundance GAMM.

Figure 2.

Figure 2. Estimated terms in GAMM model of total larval abundance, where the black line is equal to the mean effect and the grey ribbon is the 95% Confidence Interval. A) Partial linear effect of exp(-Bathymetry) in kilometres. B-D) Estimated smoothers for sampling depth (meters), latitude during Summer, latitude during Autumn, latitude during Winter, and latitude during Spring. n = 3006.

# plot Bathymetry effect
z <- evaluate_parametric_term(model2i, "exp(-Bathy)")
#z


pbath <- ggplot(z, aes(x = value, y = partial)) +
  geom_line() +
  geom_ribbon(data=z,aes(ymin=upper,ymax=lower),alpha=0.3) +
  xlab("exp(-Bathy)") +
  ylab("Partial effect of exp(-Bathy)")  +
  geom_rug(sides ="b", color="black") +
  theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + theme(plot.margin = margin(30,20,32,24))


#plot smooth effects and save as grob file
p1 <- as.grob(function() plot(model2i, shade = TRUE, select = 1))
#p1

p2 <- as.grob(function() plot(model2i, shade = TRUE, select = 2))
#p2

p3 <- as.grob(function() plot(model2i, shade = TRUE, select = 3))
#p3

p4 <- as.grob(function() plot(model2i, shade = TRUE, select = 4))
#p4

p5 <- as.grob(function() plot(model2i, shade = TRUE, select = 5))
#p5

figure_2 <- plot_grid(pbath, p1, p4, p2, p5, p3, labels = c('A', 'B', 'C', 'D', 'E', 'F'), label_size = 12, nrow = 3, ncol = 2)

figure_2

3) Richness GAMM

The following section contains code for the analysis of larval fish taxa richness using generalised additive mixed models (GAMM).

Model selection

# step 1) Full model
model1 <-  gam(richness ~  Season + Period +  exp(-Bathy)  +  s(stand_depth, bs =  "cs", k = 5) + s(dists_km, bs =  "cs", k = 5) + s(HadISST, bs =  "cs", k = 5)  + s(Latitude, bs =  "fs", xt = "cs", k = 5, by = Season, id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re"), data = eastcoast, family = nb(theta = NULL, link = "log"), method = "ML", fit=TRUE,  Select=TRUE)


#concurvity(model1, full=T) #checking concurvity

# step 2)  s(dists_km, bs =  "cs", k = 5) + s(HadISST, bs =  "cs", k = 5)  removed due to concurvity
model2 <-  gam(richness ~  Season + Period +  exp(-Bathy)  +  s(stand_depth, bs =  "cs", k = 5)  + s(Latitude, bs =  "fs", xt = "cs", k = 5, by = Season, id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re"), data = eastcoast, family = nb(theta = NULL, link = "log"), method = "ML", fit=TRUE,  Select=TRUE)


# Step 3) Period removed based on AIC score, checked with chi-sq test
model2i <-  gam(richness ~  Season  +  exp(-Bathy)  +  s(stand_depth, bs =  "cs", k = 5)  + s(Latitude, bs =  "fs", xt = "cs", k = 5, by = Season, id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re"), data = eastcoast, family = nb(theta = NULL, link = "log"), method = "ML", fit=TRUE,  Select=TRUE)

AIC(model2,model2i) ## AIC reduced with period removed
##               df      AIC
## model2  30.49724 19565.42
## model2i 29.82409 19561.91
#anova(model2,model2i, test = "Chisq") #non sig therefore can remove period (adds nothing to model)



# Step 4) linear Season term removed, but increases AIC and is sig dif in chi-sq test there retain Season in model
model3 <-  gam(richness ~  exp(-Bathy)  +  s(stand_depth, bs =  "cs", k = 5)  + s(Latitude, bs =  "fs", xt = "cs", k = 5, by = Season, id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re"), data = eastcoast, family = nb(theta = NULL, link = "log"), method = "ML", fit=TRUE,  Select=TRUE)


AIC(model3,model2i) ## AIC increased when season removed. ## model2i best model here. 
##               df      AIC
## model3  30.68741 19623.27
## model2i 29.82409 19561.91
anova(model3,model2i, test = "Chisq")  #sig dif therefore don't remove season
## Analysis of Deviance Table
## 
## Model 1: richness ~ exp(-Bathy) + s(stand_depth, bs = "cs", k = 5) + s(Latitude, 
##     bs = "fs", xt = "cs", k = 5, by = Season, id = 1) + offset(log(Volume_m3)) + 
##     s(Project_name, bs = "re")
## Model 2: richness ~ Season + exp(-Bathy) + s(stand_depth, bs = "cs", k = 5) + 
##     s(Latitude, bs = "fs", xt = "cs", k = 5, by = Season, id = 1) + 
##     offset(log(Volume_m3)) + s(Project_name, bs = "re")
##   Resid. Df Resid. Dev       Df Deviance Pr(>Chi)
## 1    2968.4     3360.4                           
## 2    2969.1     3358.8 -0.66848   1.5397
## therefore best model for richness appears to be model2i, removing Period based on AIC and s(dists_km, bs =  "cs", k = 5) + s(HadISST, bs =  "cs", k = 5) based on concurvity.

# Run with method = REML for estimates
model2i <-  gam(richness ~  Season  +  exp(-Bathy)  +  s(stand_depth, bs =  "cs", k = 5)  + s(Latitude, bs =  "fs", xt = "cs", k = 5, by = Season, id = 1) + offset(log(Volume_m3)) + s(Project_name, bs = "re"), data = eastcoast, family = nb(theta = NULL, link = "log"), method = "REML", fit=TRUE,  Select=TRUE)


#Best model summary
summary(model2i) 
## 
## Family: Negative Binomial(2.032) 
## Link function: log 
## 
## Formula:
## richness ~ Season + exp(-Bathy) + s(stand_depth, bs = "cs", k = 5) + 
##     s(Latitude, bs = "fs", xt = "cs", k = 5, by = Season, id = 1) + 
##     offset(log(Volume_m3)) + s(Project_name, bs = "re")
## 
## Parametric coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -3.98100    0.20233 -19.676  < 2e-16 ***
## Seasonspring -0.45834    0.07145  -6.415 1.41e-10 ***
## Seasonsummer  0.11461    0.04220   2.716  0.00661 ** 
## Seasonwinter -0.37823    0.06931  -5.457 4.85e-08 ***
## exp(-Bathy)   0.83862    0.07820  10.724  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                            edf Ref.df Chi.sq  p-value    
## s(stand_depth)           3.866      4  604.4  < 2e-16 ***
## s(Latitude):Seasonautumn 2.573      4  483.0  < 2e-16 ***
## s(Latitude):Seasonspring 2.169      4  749.8 7.61e-09 ***
## s(Latitude):Seasonsummer 2.714      4  178.0   0.0128 *  
## s(Latitude):Seasonwinter 2.233      4  435.8 1.76e-13 ***
## s(Project_name)          9.557     10  570.9  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  -0.0346   Deviance explained = 45.4%
## -REML = 9823.8  Scale est. = 1         n = 3000
#par(mfcol=c(3,2))
#plot(model2i, shade = TRUE) ### spline plots

#checking concurvity
#concurvity(model2i,full=TRUE)

#checking assumptions
appraise(model2i)

#overdisperison test
E1 <- resid(model2i, type = "pearson")   #ok - less than 2
N <- nrow(eastcoast)
p <- length(coef(model2i))
Overdispersion <- sum(E1^2) / (N - p)
#paste('Overdispersion: ',Overdispersion) #"Overdispersion:  1.01888080182827"

Checking spatial auto correlation of Taxa Richness GAMM with Moran’s I If Moran’s I is ~0 then there is no spatial autocorrelation.

res <- residuals(model2i)
#res


Coords <-cbind(eastcoast$Latitude, eastcoast$Longitude)

bw <- 6

mI <-moransI(Coords,bw,res)

moran.table <-matrix(data=NA,nrow=1,ncol=6)

col.names <-c("Moran's I", "Expected I", "Z resampling", "P-value resampling","Z randomization", "P-value randomization")

colnames(moran.table) <- col.names

moran.table[1,1] <- mI$Morans.I

moran.table[1,2] <- mI$Expected.I

moran.table[1,3] <- mI$z.resampling

moran.table[1,4] <- mI$p.value.resampling

moran.table[1,5] <- mI$z.randomization

moran.table[1,6] <- mI$p.value.randomization

print(moran.table)

Plotting output

Plotting effects of the larval fish taxa richness GAMM.

Figure 3.

Figure 3. Estimated terms in GAMM model of larval taxa richness, where the black line is equal to the mean effect and the grey ribbon is 95% Confidence Interval. A) Partial linear effect of exp(-Bathymetry) in kilometres. B) Non-linear effect of sampling depth (meters). C) Non-linear effect of latitude during Summer. D) Non-linear effect of latitude during Autumn. E) Non-linear effect of latitude during Winter. F) Non-linear effect of latitude during Spring. n = 3006.

z <- evaluate_parametric_term(model2i, "exp(-Bathy)")
#z


pbath <- ggplot(z, aes(x = value, y = partial)) +
  geom_line() +
  geom_ribbon(data=z,aes(ymin=upper,ymax=lower),alpha=0.3) +
  xlab("exp(-Bathy)") +
  ylab("Partial effect of exp(-Bathy)")  +
  geom_rug(sides ="b", color="black") +
  theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + theme(plot.margin = margin(30,20,32,24))
  


p1 <- as.grob(function() plot(model2i, shade = TRUE, select = 1))
#p1

p2 <- as.grob(function() plot(model2i, shade = TRUE, select = 2))
#p2


p3 <- as.grob(function() plot(model2i, shade = TRUE, select = 3))
#p3

p4 <- as.grob(function() plot(model2i, shade = TRUE, select = 4))
#p4

p5 <- as.grob(function() plot(model2i, shade = TRUE, select = 5))
#p5

figure_3 <- plot_grid(pbath, p1, p4, p2, p5, p3, labels = c('A', 'B', 'C', 'D', 'E', 'F'), label_size = 12, nrow = 3, ncol = 2)

figure_3

4) Multivariate Analyses

This section of the scripts uses an intercept-only multivariate generalised linear model (MGLM) model of the larval fish community in a Gaussian copula graphical model (GCGM), and quantile regression of the , to observe the effects of temperature (Figure 4), latitude (Supplementary Figure 2) and time period (Supplementary Figure 3).

Select inshore species in matrix to analyse

elements <- inshore[1:2428,22:239]



myvars <- names(elements) %in% c("Acropomatidae_Synagrops.spp_37311949", 
                            "Aploactinidae_other_37290000",
                            "Bothidae_Crossorhombus.spp_37460907",
                            "Bovichtidae_Pseudaphritis.urvilli_37403003",
                            "Callanthiidae_other_37311957",
                            "Cepolidae_Acanthocepola.spp_37380901",
                            "Cepolidae_Owstonia.spp_37380903",
                            "Cepolidae_other_37380000",
                            "Cetomimidae_37132000",
                            "Chandidae_other_37310900",
                            "Engraulidae_other_37086000",
                            "Ipnopidae_37123000",
                            "Latridae_other_37378000",
                            "Leptobramidae_37357905",
                            "Microcanthidae_other_37361900",
                            "Ophidiidae_Brotula.spp_37228912",
                            "Plesiopidae_37316000",
                            "Terapontidae_other_37321000",
                            "Xiphiidae_Xiphias.gladius_37442001") 


elements <- elements[!myvars] #removing species which never occur on in less than 200m bathymetry on the east coast.


elements <- mvabund(elements[,])

Fitting GCGM with intercept only MGLM

#boxplot(elements) # inspect data ranges

#meanvar.plot(elements) # check mean variance relationship


#fit <- manyglm(elements ~Latitude, family = "negative.binomial", offset = log(Volume_m3),data=inshore)

fit_null <- manyglm(elements ~1, family = "negative.binomial", offset = log(Volume_m3),data=inshore)



fit_cop <- cord(fit_null)
plot(fit_cop,biplot=TRUE)

score_1 <- fit_cop$scores[[1]][1,]
score_2 <- fit_cop$scores[[1]][2,]

Plot GCGM bivariate scores 1 and 2 with SST, and then latitude.

plot_dat <- data.frame(score_1,score_2,inshore$Region, inshore$Latitude, inshore$HadISST, inshore$Period)

SST_bivariate_plot <- ggplot(plot_dat,aes(x=score_1,y=score_2,col=inshore.HadISST))+
  geom_point() +
  xlab("Score 1") +
  ylab("Score 2") +
  scale_color_continuous(high = "red", low = "blue") +
  theme_bw() +
  theme(legend.position=c(0.85, .9), legend.title = element_text(size=10, face="bold"),
               legend.text = element_text(size = 10, face = "bold"))
  


LAT_bivariate_plot <- ggplot(plot_dat,aes(x=score_1,y=score_2,col=inshore.Latitude))+
  geom_point() +
  xlab("Score 1") +
  ylab("Score 2") +
  scale_color_continuous(high = "red", low = "blue") +
  theme_bw() +
  theme(legend.position=c(0.85, .9), legend.title = element_text(size=10, face="bold"),
               legend.text = element_text(size = 10, face = "bold"))

Quantile regression of bivariate scores against SST

Quantile regression of score 1 from GCGM against SST

qs <- c(0.05, 0.5, 0.95) ### What quantiles are we interested in... 5%, 50% and 95%

#quantile regression for score 1 
qr1 <- rq(score_1 ~ inshore.HadISST, data=plot_dat, tau = qs)

SM <- summary(qr1, se = "boot", bsmethod = "wild", R = 1000)  ## method of calculating SE and bootstrap method if applicable
SM
## 
## Call: rq(formula = score_1 ~ inshore.HadISST, tau = qs, data = plot_dat)
## 
## tau: [1] 0.05
## 
## Coefficients:
##                 Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      -1.46418   0.08469  -17.28774   0.00000
## inshore.HadISST   0.02291   0.00444    5.15674   0.00000
## 
## Call: rq(formula = score_1 ~ inshore.HadISST, tau = qs, data = plot_dat)
## 
## tau: [1] 0.5
## 
## Coefficients:
##                 Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      -1.50302   0.09426  -15.94554   0.00000
## inshore.HadISST   0.06436   0.00526   12.23176   0.00000
## 
## Call: rq(formula = score_1 ~ inshore.HadISST, tau = qs, data = plot_dat)
## 
## tau: [1] 0.95
## 
## Coefficients:
##                 Value    Std. Error t value  Pr(>|t|)
## (Intercept)     -4.01042  0.53902   -7.44027  0.00000
## inshore.HadISST  0.29237  0.02887   10.12823  0.00000

Plotting quantile regression score 1 from GCGM against SST

### Build a ggplot
p <- ggplot(plot_dat, aes(x = inshore.HadISST, y = score_1)) +
  geom_point(size = 3, na.rm = TRUE, shape = 1) 
   
p <- p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

p <- p + geom_abline(intercept = qr1$coefficients[1,1], slope = qr1$coefficients[2,1], 
                     linetype = 'longdash', col = "blue") # 5th quantile
p <- p + geom_abline(intercept = qr1$coefficients[1,2], slope = qr1$coefficients[2,2], col = "blue")
                   #  linetype = 'longdash' ) # 50th quantile
p <- p + geom_abline(intercept = qr1$coefficients[1,3], slope = qr1$coefficients[2,3], 
                     linetype = 'longdash', col = "blue") # 95th quantile
p <- p + xlab("SST (°C)")
p <- p + ylab("Score 1")
p <- p + theme(axis.text=element_text(size=16, face = "bold"))
p <- p + theme(axis.title=element_text(size=16, face = "bold"))
p_score1 <- p + theme(legend.position=c(0.15, .8), legend.title = element_text(size=16, face="bold"),
               legend.text = element_text(size = 16, face = "bold"))
#p_score1

Quantile regression of score 2 from GCGM against SST

qs <- c(0.05, 0.5, 0.95) ### What quantiles are we interested in... 5%, 50% and 95%

#quantile regression for score 1 
qr2 <- rq(score_2 ~ inshore.HadISST, data=plot_dat, tau = qs)

SM <- summary(qr2, se = "boot", bsmethod = "wild", R = 1000)  ## method of calculating SE and bootstrap method if applicable
SM
## 
## Call: rq(formula = score_2 ~ inshore.HadISST, tau = qs, data = plot_dat)
## 
## tau: [1] 0.05
## 
## Coefficients:
##                 Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      -1.84492   0.13320  -13.85082   0.00000
## inshore.HadISST   0.04684   0.00723    6.47927   0.00000
## 
## Call: rq(formula = score_2 ~ inshore.HadISST, tau = qs, data = plot_dat)
## 
## tau: [1] 0.5
## 
## Coefficients:
##                 Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      -2.17395   0.12361  -17.58653   0.00000
## inshore.HadISST   0.10387   0.00657   15.80505   0.00000
## 
## Call: rq(formula = score_2 ~ inshore.HadISST, tau = qs, data = plot_dat)
## 
## tau: [1] 0.95
## 
## Coefficients:
##                 Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      -4.28122   0.13235  -32.34750   0.00000
## inshore.HadISST   0.28888   0.00879   32.86453   0.00000

Plotting quantile regression score 2 from GCGM against SST

### Build a ggplot
p <- ggplot(plot_dat, aes(x = inshore.HadISST, y = score_2)) +
  geom_point(size = 3, na.rm = TRUE, shape = 1) 
p <- p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) ## set up plot area
p <- p + geom_abline(intercept = qr2$coefficients[1,1], slope = qr2$coefficients[2,1], 
                     linetype = 'longdash', col = "blue") # 5th quantile
p <- p + geom_abline(intercept = qr2$coefficients[1,2], slope = qr2$coefficients[2,2], col = "blue")
                   #  linetype = 'longdash' ) # 50th quantile
p <- p + geom_abline(intercept = qr2$coefficients[1,3], slope = qr2$coefficients[2,3], 
                     linetype = 'longdash', col = "blue") # 95th quantile
p <- p + xlab("SST (°C)")
p <- p + ylab("Score 2")
p <- p + theme(axis.text=element_text(size=16, face = "bold"))
p <- p + theme(axis.title=element_text(size=16, face = "bold"))
p_score2 <- p + theme(legend.position=c(0.15, .8), legend.title = element_text(size=16, face="bold"),
               legend.text = element_text(size = 16, face = "bold"))
#p_score2

Plotting the GCGM bivariate scores 1 and 2 with SST, and quantile regression plots, as per Figure 4 in the manuscript.

Figure 4.

Figure 4. A) GCGM bivariate scores 1 and 2 from the MGLM intercept-only model of all larval fish taxa abundances . Each point represents an assemblage and the colour of the points represents SST in °C. B) GCGM bivariate score 1 by SST and C) GCGM bivariate score 2 by SST. Blue lines represent 0.05, 0.5 and 0.95 quantiles. n = 2428.

figure_4 <- plot_grid(SST_bivariate_plot, p_score1, p_score2, labels = c('A', 'B', 'C'), label_size = 12, nrow = 1, ncol = 3)

figure_4

Univariate tests for MGLM with SST

The univariate tests for the MGLM with SST were run on the University of New South Wales computational cluster ‘Katana’ supported by Research Technology Services. The r script included in the repo ‘manyglm_uni_test.R’ contains the analyses. Here, we import the output to observe the results.

Read in rds and display output

Uni_test_manyglm_output <- readRDS("outputs/Uni_test_manyglm_output.rds")

#Uni_test_manyglm_output

#Uni_test_manyglm_output$uni.p
uni_p <- as.data.frame(Uni_test_manyglm_output$uni.p)
#uni_p


uni_ts <- as.data.frame(Uni_test_manyglm_output$uni.test)
#uni_ts

useful_output <- rbind(uni_p, uni_ts)
#useful_output

useful_output <- t(useful_output)

useful_output <- as.data.frame(useful_output)

sig_out <- subset(useful_output, useful_output$HadISST <= 0.05)

#sig_out




#sum(useful_output$HadISST1) #whole test dev

useful_output$dev_exp <- (useful_output$HadISST1/sum(useful_output$HadISST1))*100

#useful_output


useful_output <- useful_output %>% 
  dplyr::rename(
    P = HadISST,
    TS = HadISST1
    ) %>%
  dplyr::select(P, TS, dev_exp)



head(useful_output)
##                                                P         TS      dev_exp
## Acanthuridae_37437000                     0.0001 156.978730 1.1467679622
## Acropomatidae_Apogonops.anomalus_37311053 0.0001  95.979510 0.7011537639
## Acropomatidae_other_37311956              1.0000   0.111148 0.0008119635
## Ammodytidae_37425000                      0.0002  56.514417 0.4128516210
## Anguilliformes_37990019                   0.0001 271.926545 1.9864898281
## Antennariidae_37210000                    0.4399   9.786945 0.0714960249
#write.csv(useful_output, "manyglm_temp_uni.csv")

Quantile regression of bivariate scores against Latitude

Now running quantile regression of score 1 from GCGM against Latitude

qs <- c(0.05, 0.5, 0.95) ### What quantiles are we interested in... 5%, 50% and 95%

#quantile regression for score 1 
qr1 <- rq(score_1 ~ inshore.Latitude, data=plot_dat, tau = qs)

SM <- summary(qr1, se = "boot", bsmethod = "wild", R = 1000)  ## method of calculating SE and bootstrap method if applicable
SM
## 
## Call: rq(formula = score_1 ~ inshore.Latitude, tau = qs, data = plot_dat)
## 
## tau: [1] 0.05
## 
## Coefficients:
##                  Value    Std. Error t value  Pr(>|t|)
## (Intercept)      -0.46350  0.17303   -2.67868  0.00744
## inshore.Latitude  0.01580  0.00468    3.37658  0.00075
## 
## Call: rq(formula = score_1 ~ inshore.Latitude, tau = qs, data = plot_dat)
## 
## tau: [1] 0.5
## 
## Coefficients:
##                  Value    Std. Error t value  Pr(>|t|)
## (Intercept)       1.42931  0.11535   12.39081  0.00000
## inshore.Latitude  0.04724  0.00306   15.45606  0.00000
## 
## Call: rq(formula = score_1 ~ inshore.Latitude, tau = qs, data = plot_dat)
## 
## tau: [1] 0.95
## 
## Coefficients:
##                  Value    Std. Error t value  Pr(>|t|)
## (Intercept)       9.58603  0.40829   23.47826  0.00000
## inshore.Latitude  0.22208  0.01057   21.01611  0.00000

Plotting quantile regression score 1 from GCGM against Latitude

### Build a ggplot
p <- ggplot(plot_dat, aes(x = inshore.Latitude, y = score_1)) +
  geom_point(size = 3, na.rm = TRUE, shape = 1) 
   
p <- p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

p <- p + geom_abline(intercept = qr1$coefficients[1,1], slope = qr1$coefficients[2,1], 
                     linetype = 'longdash', col = "blue") # 5th quantile
p <- p + geom_abline(intercept = qr1$coefficients[1,2], slope = qr1$coefficients[2,2], col = "blue")
                   #  linetype = 'longdash' ) # 50th quantile
p <- p + geom_abline(intercept = qr1$coefficients[1,3], slope = qr1$coefficients[2,3], 
                     linetype = 'longdash', col = "blue") # 95th quantile
p <- p + xlab("Latitude (°)")
p <- p + ylab("Score 1")
p <- p + theme(axis.text=element_text(size=16, face = "bold"))
p <- p + theme(axis.title=element_text(size=16, face = "bold"))
p_score1 <- p + theme(legend.position=c(0.15, .8), legend.title = element_text(size=16, face="bold"),
               legend.text = element_text(size = 16, face = "bold"))
#p_score1

Now running quantile regression of score 1 from GCGM against Latitude

qs <- c(0.05, 0.5, 0.95) ### What quantiles are we interested in... 5%, 50% and 95%

#quantile regression for score 1 
qr2 <- rq(score_2 ~ inshore.Latitude, data=plot_dat, tau = qs)

SM <- summary(qr2, se = "boot", bsmethod = "wild", R = 1000)  ## method of calculating SE and bootstrap method if applicable
SM
## 
## Call: rq(formula = score_2 ~ inshore.Latitude, tau = qs, data = plot_dat)
## 
## tau: [1] 0.05
## 
## Coefficients:
##                  Value    Std. Error t value  Pr(>|t|)
## (Intercept)      -0.08068  0.21772   -0.37058  0.71098
## inshore.Latitude  0.02462  0.00601    4.09674  0.00004
## 
## Call: rq(formula = score_2 ~ inshore.Latitude, tau = qs, data = plot_dat)
## 
## tau: [1] 0.5
## 
## Coefficients:
##                  Value    Std. Error t value  Pr(>|t|)
## (Intercept)       1.58956  0.16467    9.65318  0.00000
## inshore.Latitude  0.04990  0.00438   11.38512  0.00000
## 
## Call: rq(formula = score_2 ~ inshore.Latitude, tau = qs, data = plot_dat)
## 
## tau: [1] 0.95
## 
## Coefficients:
##                  Value    Std. Error t value  Pr(>|t|)
## (Intercept)       8.61369  0.25249   34.11502  0.00000
## inshore.Latitude  0.19815  0.00614   32.26831  0.00000

Plotting quantile regression score 2 from GCGM against Latitude

### Buiild a ggplot
p <- ggplot(plot_dat, aes(x = inshore.Latitude, y = score_2)) +
  geom_point(size = 3, na.rm = TRUE, shape = 1) 
p <- p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) ## set up plot area
p <- p + geom_abline(intercept = qr2$coefficients[1,1], slope = qr2$coefficients[2,1], 
                     linetype = 'longdash', col = "blue") # 5th quantile
p <- p + geom_abline(intercept = qr2$coefficients[1,2], slope = qr2$coefficients[2,2], col = "blue")
                   #  linetype = 'longdash' ) # 50th quantile
p <- p + geom_abline(intercept = qr2$coefficients[1,3], slope = qr2$coefficients[2,3], 
                     linetype = 'longdash', col = "blue") # 95th quantile
p <- p + xlab("Latitude (°)")
p <- p + ylab("Score 2")
p <- p + theme(axis.text=element_text(size=16, face = "bold"))
p <- p + theme(axis.title=element_text(size=16, face = "bold"))
p_score2 <- p + theme(legend.position=c(0.15, .8), legend.title = element_text(size=16, face="bold"),
               legend.text = element_text(size = 16, face = "bold"))
#p_score2

Plotting the GCGM bivariate scores 1 and 2 with Latitude, and quantile regression plots, as per Supplementary Figure 2 in the manuscript.

Supp Figure 2.

Supplementary Figure 2. A) GCGM bivariate scores 1 and 2. Each point represents an assemblage and the colour of the points represents latitude. B) GCGM bivariate score 1 by latitude and C) GCGM bivariate score 2 by latitude. Blue lines represent 0.05, 0.5 and 0.95 quantiles. n = 2428.

supp_figure_2 <- plot_grid(LAT_bivariate_plot, p_score1, p_score2, labels = c('A', 'B', 'C'), label_size = 12, nrow = 1, ncol = 3)

supp_figure_2

Observing difference in community composition of larval fish taxa with latitude before and after 1998 (coinciding with the 1998/99 El Nino event)in the bivariate scores of the GCGM with intercept only MGLM.

Figure 5.

Figure 5. A) GCGM bivariate scores 1 and 2. Each point represents an assemblage and the colour of the points represents the ‘Period’, black circles are pre-1998 and orange circles are post-1998. B) GCGM bivariate score 1 by latitude and C) GCGM bivariate score 2 by latitude. Trend lines are fit with a loess smoother.

my_cols <- c("#E69F00", "#000000")

p1<- ggplot(plot_dat,aes(x=score_1,y=score_2))+
  geom_point(aes(color = inshore.Period), size = 4, alpha = 0.5) + 
  scale_color_manual(values = my_cols) +
  theme_classic() +
  theme(legend.position=c(0.85, .9), legend.title = element_text(size=10, face="bold"),
               legend.text = element_text(size = 10, face = "bold"))
  
#p1

p2 <- ggplot(plot_dat,aes(x=inshore.Latitude, y=score_1, col=inshore.Period))+
   geom_point(aes(color = inshore.Period), size = 4, alpha = 0.5,  shape = 1) + 
  scale_color_manual(values = my_cols) +
  theme_classic() +
  xlab("Latitude (°C)") +
  theme(legend.position="none") +
  geom_smooth(method = "loess")


p3 <-  ggplot(plot_dat,aes(x=inshore.Latitude, y=score_2, col=inshore.Period))+
  geom_point(aes(color = inshore.Period), size = 4, alpha = 0.5,  shape = 1) + 
  scale_color_manual(values = my_cols) +
  theme_classic() +
  xlab("Latitude (°C)") +
  theme(legend.position="none") +
  geom_smooth(method = "loess")



supp_figure_3 <- plot_grid(p1, p2, p3, labels = c('A', 'B', 'C'), label_size = 12, nrow = 1, ncol = 3)

supp_figure_3

Ssssion Info

devtools::session_info()
## ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 3.6.1 (2019-07-05)
##  os       macOS Mojave 10.14.6        
##  system   x86_64, darwin15.6.0        
##  ui       X11                         
##  language (EN)                        
##  collate  en_AU.UTF-8                 
##  ctype    en_AU.UTF-8                 
##  tz       Australia/Sydney            
##  date     2020-07-07                  
## 
## ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
##  package      * version    date       lib source                              
##  abind          1.4-5      2016-07-21 [1] CRAN (R 3.6.0)                      
##  acepack        1.4.1      2016-10-29 [1] CRAN (R 3.6.0)                      
##  assertthat     0.2.1      2019-03-21 [1] CRAN (R 3.6.0)                      
##  backports      1.1.5      2019-10-02 [1] CRAN (R 3.6.0)                      
##  base64enc      0.1-3      2015-07-28 [1] CRAN (R 3.6.0)                      
##  BiocManager    1.30.9     2019-10-23 [1] CRAN (R 3.6.1)                      
##  broom          0.5.2      2019-04-07 [1] CRAN (R 3.6.0)                      
##  callr          3.3.2      2019-09-22 [1] CRAN (R 3.6.0)                      
##  car          * 3.0-5      2019-11-15 [1] CRAN (R 3.6.0)                      
##  carData      * 3.0-3      2019-11-16 [1] CRAN (R 3.6.0)                      
##  cellranger     1.1.0      2016-07-27 [1] CRAN (R 3.6.0)                      
##  checkmate      1.9.4      2019-07-04 [1] CRAN (R 3.6.0)                      
##  cli            2.0.2      2020-02-28 [1] CRAN (R 3.6.0)                      
##  cluster        2.1.0      2019-06-19 [1] CRAN (R 3.6.1)                      
##  colorspace     1.4-1      2019-03-18 [1] CRAN (R 3.6.0)                      
##  cowplot      * 1.0.0      2019-07-11 [1] CRAN (R 3.6.0)                      
##  crayon         1.3.4      2017-09-16 [1] CRAN (R 3.6.0)                      
##  curl           4.3        2019-12-02 [1] CRAN (R 3.6.0)                      
##  data.table     1.12.8     2019-12-09 [1] CRAN (R 3.6.0)                      
##  DBI            1.1.0      2019-12-15 [1] CRAN (R 3.6.0)                      
##  dbplyr         1.4.2      2019-06-17 [1] CRAN (R 3.6.0)                      
##  desc           1.2.0      2018-05-01 [1] CRAN (R 3.6.0)                      
##  devtools       2.2.1      2019-09-24 [1] CRAN (R 3.6.0)                      
##  digest         0.6.25     2020-02-23 [1] CRAN (R 3.6.0)                      
##  dplyr        * 1.0.0      2020-05-29 [1] CRAN (R 3.6.2)                      
##  ecoCopula    * 0.0.0.9000 2020-07-03 [1] Github (gordy2x/ecoCopula@ea048e1)  
##  ellipsis       0.3.0      2019-09-20 [1] CRAN (R 3.6.0)                      
##  evaluate       0.14       2019-05-28 [1] CRAN (R 3.6.0)                      
##  fansi          0.4.1      2020-01-08 [1] CRAN (R 3.6.0)                      
##  farver         2.0.1      2019-11-13 [1] CRAN (R 3.6.0)                      
##  forcats      * 0.4.0      2019-02-17 [1] CRAN (R 3.6.0)                      
##  foreign        0.8-72     2019-08-02 [1] CRAN (R 3.6.0)                      
##  Formula        1.2-3      2018-05-03 [1] CRAN (R 3.6.0)                      
##  fs             1.3.1      2019-05-06 [1] CRAN (R 3.6.0)                      
##  gdata          2.18.0     2017-06-06 [1] CRAN (R 3.6.0)                      
##  generics       0.0.2      2018-11-29 [1] CRAN (R 3.6.0)                      
##  ggplot2      * 3.2.1      2019-08-10 [1] CRAN (R 3.6.0)                      
##  ggplotify    * 0.0.5      2020-03-12 [1] CRAN (R 3.6.0)                      
##  ggpubr       * 0.2.3      2019-09-03 [1] CRAN (R 3.6.0)                      
##  ggsignif       0.6.0      2019-08-08 [1] CRAN (R 3.6.0)                      
##  glue           1.4.0      2020-04-03 [1] CRAN (R 3.6.1)                      
##  gratia       * 0.2-3      2020-07-01 [1] Github (gavinsimpson/gratia@245ac45)
##  gridExtra      2.3        2017-09-09 [1] CRAN (R 3.6.0)                      
##  gridGraphics * 0.5-0      2020-02-25 [1] CRAN (R 3.6.0)                      
##  gtable         0.3.0      2019-03-25 [1] CRAN (R 3.6.0)                      
##  gtools         3.8.1      2018-06-26 [1] CRAN (R 3.6.0)                      
##  haven          2.2.0      2019-11-08 [1] CRAN (R 3.6.0)                      
##  Hmisc          4.4-0      2020-03-23 [1] CRAN (R 3.6.0)                      
##  hms            0.5.3      2020-01-08 [1] CRAN (R 3.6.0)                      
##  htmlTable      1.13.3     2019-12-04 [1] CRAN (R 3.6.0)                      
##  htmltools      0.4.0      2019-10-04 [1] CRAN (R 3.6.0)                      
##  htmlwidgets    1.5.1      2019-10-08 [1] CRAN (R 3.6.0)                      
##  httr           1.4.1      2019-08-05 [1] CRAN (R 3.6.0)                      
##  jpeg           0.1-8      2014-01-23 [1] CRAN (R 3.6.0)                      
##  jsonlite       1.6.1      2020-02-02 [1] CRAN (R 3.6.0)                      
##  knitr          1.26       2019-11-12 [1] CRAN (R 3.6.0)                      
##  labeling       0.3        2014-08-23 [1] CRAN (R 3.6.0)                      
##  lattice      * 0.20-38    2018-11-04 [1] CRAN (R 3.6.1)                      
##  latticeExtra   0.6-29     2019-12-19 [1] CRAN (R 3.6.0)                      
##  lazyeval       0.2.2      2019-03-15 [1] CRAN (R 3.6.0)                      
##  lctools      * 0.2-8      2020-04-06 [1] CRAN (R 3.6.2)                      
##  lifecycle      0.2.0      2020-03-06 [1] CRAN (R 3.6.0)                      
##  lubridate    * 1.7.4      2018-04-11 [1] CRAN (R 3.6.0)                      
##  magrittr     * 1.5        2014-11-22 [1] CRAN (R 3.6.0)                      
##  MASS           7.3-51.4   2019-03-31 [1] CRAN (R 3.6.1)                      
##  Matrix         1.2-17     2019-03-22 [1] CRAN (R 3.6.1)                      
##  MatrixModels   0.4-1      2015-08-22 [1] CRAN (R 3.6.0)                      
##  memoise        1.1.0      2017-04-21 [1] CRAN (R 3.6.0)                      
##  mgcv         * 1.8-28     2019-03-21 [1] CRAN (R 3.6.1)                      
##  mice           3.8.0      2020-02-21 [1] CRAN (R 3.6.0)                      
##  modelr         0.1.5      2019-08-08 [1] CRAN (R 3.6.0)                      
##  munsell        0.5.0      2018-06-12 [1] CRAN (R 3.6.0)                      
##  mvabund      * 4.1.3      2020-02-27 [1] CRAN (R 3.6.0)                      
##  mvtnorm        1.0-11     2019-06-19 [1] CRAN (R 3.6.0)                      
##  nlme         * 3.1-141    2019-08-01 [1] CRAN (R 3.6.0)                      
##  nnet           7.3-12     2016-02-02 [1] CRAN (R 3.6.1)                      
##  openxlsx       4.1.4      2019-12-06 [1] CRAN (R 3.6.0)                      
##  permute      * 0.9-5      2019-03-12 [1] CRAN (R 3.6.0)                      
##  pillar         1.4.3      2019-12-20 [1] CRAN (R 3.6.0)                      
##  pkgbuild       1.0.6      2019-10-09 [1] CRAN (R 3.6.0)                      
##  pkgconfig      2.0.3      2019-09-22 [1] CRAN (R 3.6.0)                      
##  pkgload        1.0.2      2018-10-29 [1] CRAN (R 3.6.0)                      
##  plyr           1.8.6      2020-03-03 [1] CRAN (R 3.6.0)                      
##  png            0.1-7      2013-12-03 [1] CRAN (R 3.6.0)                      
##  pracma       * 2.2.5      2019-04-09 [1] CRAN (R 3.6.0)                      
##  prettyunits    1.1.1      2020-01-24 [1] CRAN (R 3.6.0)                      
##  processx       3.4.1      2019-07-18 [1] CRAN (R 3.6.0)                      
##  ps             1.3.0      2018-12-21 [1] CRAN (R 3.6.0)                      
##  pscl           1.5.5      2020-03-07 [1] CRAN (R 3.6.0)                      
##  purrr        * 0.3.4      2020-04-17 [1] CRAN (R 3.6.2)                      
##  quantreg     * 5.52       2019-11-09 [1] CRAN (R 3.6.0)                      
##  R6             2.4.1      2019-11-12 [1] CRAN (R 3.6.0)                      
##  RColorBrewer   1.1-2      2014-12-07 [1] CRAN (R 3.6.0)                      
##  Rcpp           1.0.4.6    2020-04-09 [1] CRAN (R 3.6.1)                      
##  readr        * 1.3.1      2018-12-21 [1] CRAN (R 3.6.0)                      
##  readxl         1.3.1      2019-03-13 [1] CRAN (R 3.6.0)                      
##  remotes        2.1.0      2019-06-24 [1] CRAN (R 3.6.0)                      
##  reprex         0.3.0      2019-05-16 [1] CRAN (R 3.6.0)                      
##  reshape        0.8.8      2018-10-23 [1] CRAN (R 3.6.0)                      
##  rio            0.5.16     2018-11-26 [1] CRAN (R 3.6.0)                      
##  rlang          0.4.6      2020-05-02 [1] CRAN (R 3.6.2)                      
##  rmarkdown      2.1        2020-01-20 [1] CRAN (R 3.6.0)                      
##  rpart          4.1-15     2019-04-12 [1] CRAN (R 3.6.1)                      
##  rprojroot      1.3-2      2018-01-03 [1] CRAN (R 3.6.0)                      
##  rstudioapi     0.10       2019-03-19 [1] CRAN (R 3.6.0)                      
##  rvcheck        0.1.8      2020-03-01 [1] CRAN (R 3.6.0)                      
##  rvest          0.3.5      2019-11-08 [1] CRAN (R 3.6.0)                      
##  scales         1.1.0      2019-11-18 [1] CRAN (R 3.6.0)                      
##  sessioninfo    1.1.1      2018-11-05 [1] CRAN (R 3.6.0)                      
##  sp             1.3-2      2019-11-07 [1] CRAN (R 3.6.0)                      
##  SparseM      * 1.77       2017-04-23 [1] CRAN (R 3.6.0)                      
##  statmod        1.4.34     2020-02-17 [1] CRAN (R 3.6.0)                      
##  stringi        1.4.6      2020-02-17 [1] CRAN (R 3.6.0)                      
##  stringr      * 1.4.0      2019-02-10 [1] CRAN (R 3.6.0)                      
##  survival       3.1-12     2020-04-10 [1] CRAN (R 3.6.2)                      
##  testthat       2.2.1      2019-07-25 [1] CRAN (R 3.6.0)                      
##  tibble       * 3.0.1      2020-04-20 [1] CRAN (R 3.6.2)                      
##  tidyr        * 1.1.0      2020-05-20 [1] CRAN (R 3.6.2)                      
##  tidyselect     1.1.0      2020-05-11 [1] CRAN (R 3.6.2)                      
##  tidyverse    * 1.3.0      2019-11-21 [1] CRAN (R 3.6.0)                      
##  tweedie        2.3.2      2017-12-14 [1] CRAN (R 3.6.0)                      
##  usethis        1.5.1      2019-07-04 [1] CRAN (R 3.6.0)                      
##  vctrs          0.3.1      2020-06-05 [1] CRAN (R 3.6.2)                      
##  vegan        * 2.5-6      2019-09-01 [1] CRAN (R 3.6.0)                      
##  weights        1.0.1      2020-02-12 [1] CRAN (R 3.6.0)                      
##  withr          2.1.2      2018-03-15 [1] CRAN (R 3.6.0)                      
##  xfun           0.11       2019-11-12 [1] CRAN (R 3.6.0)                      
##  xml2           1.3.0      2020-04-01 [1] CRAN (R 3.6.2)                      
##  yaml           2.2.1      2020-02-01 [1] CRAN (R 3.6.0)                      
##  zip            2.0.4      2019-09-01 [1] CRAN (R 3.6.0)                      
## 
## [1] /Library/Frameworks/R.framework/Versions/3.6/Resources/library